/* "Efficiency and water use: Dynamic effects of irrigation technology adoption"
by Micah Cameron-Harp and Nathan Hendricks

Code written by Micah Cameron-Harp
May 16th, 2024

This do file creates Figures A4 and A5 displaying the pre-treatment irrigation
	behavior by cohort.
*/

*Define directories and set working directory
/* NOTE - To replicate our results, you need to change the root directory address in the next line */
global dr_root = "C:\Users\Micah\Dropbox\Irrigation technology transition\final revisions for conditional acceptance\replication materials"
global dr_code = "${dr_root}\code"
global dr_data = "${dr_root}\data"
global dr_output = "${dr_root}\outputs"
global dr_output_main = "${dr_root}\outputs\main_text"
global dr_output_app = "${dr_root}\outputs\appendices"
global dr_output_log = "${dr_root}\outputs\logs"
global dr_temp = "${dr_root}\data\intermediate"
cd "${dr_root}"

*********************** Figure A.4 *********************************************
*Load in the water right group collapsed data for the flood to cp/lepa transition
import delimited using "$dr_temp\flood_cporlepa_prepped.csv", clear
*Limit to period of analysis
keep if inrange(wua_year, 1991, 2015) 
*Make mean of each variables
foreach var of varlist af_used acres_irr depth_applied {
egen mean_`var' = mean(`var')
}
*Drop not-yet-treated during the times when they are treated
drop if wua_year >= cohort_flood_cplepa & cohort_flood_cplepa!=0

*Create long difference variables
xtset wr_group wua_year

*Keep minimum year in dataset and last year before treatment
keep if wua_year==min_year | wua_year==cohort_flood_cplepa - 1
drop if min_year==cohort_flood_cplepa - 1
gen ind_min_year = 0
replace ind_min_year = 1 if min_year==wua_year

foreach var of varlist af_used acres_irr depth_applied {
	gen min_year_`var' = `var'*ind_min_year
	bysort wr_group: egen min_`var' = max(min_year_`var')
	drop min_year_`var'
}
keep if wua_year==cohort_flood_cplepa - 1

foreach var of varlist af_used acres_irr depth_applied {
	gen diff_`var' = `var' - min_`var'
	replace diff_`var' = 100*diff_`var'/mean_`var'
}

foreach var of varlist diff_af_used diff_acres_irr diff_depth_applied {
* Use egen to generate the median, quartiles, interquartile range (IQR), and mean. 
bysort cohort_flood_cplepa: egen med_`var' = median(`var')
by cohort_flood_cplepa: egen lqt_`var' = pctile(`var'), p(25)
by cohort_flood_cplepa: egen uqt_`var' = pctile(`var'), p(75)
by cohort_flood_cplepa: egen iqr_`var' = iqr(`var')
by cohort_flood_cplepa: egen mean_`var' = mean(`var')

* Find the lowest value that is more than lqt - 1.5 iqr 
* this is used to form the lower "whisker" of the boxplot.
gen l_`var' = `var' if(`var' >= lqt_`var'-1.5*iqr_`var')
by cohort_flood_cplepa: egen ls_`var' = min(l_`var')

* Find the highest value that is less than uqt + 1.5 iqr 
* this is used to form the upper "whisker" of the boxplot.
gen u_`var' = `var' if(`var' <= uqt_`var'+1.5*iqr_`var')
by cohort_flood_cplepa: egen us_`var' = max(u_`var')

* Find any outliers (i.e. values more than 1.5 IQRs from the upper and lower quartiles). 
gen outliers_`var' = `var' if(`var'<=lqt_`var'-1.5*iqr_`var' | `var'>=uqt_`var'+1.5*iqr_`var')
}


*Plot for af_used
twoway rbar lqt_diff_af_used med_diff_af_used cohort_flood_cplepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rbar med_diff_af_used uqt_diff_af_used cohort_flood_cplepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rspike lqt_diff_af_used ls_diff_af_used cohort_flood_cplepa, lwidth(.1) lcolor(black) || ///
       rspike uqt_diff_af_used us_diff_af_used cohort_flood_cplepa, lwidth(.1) lcolor(black) || ///
       rcap ls_diff_af_used ls_diff_af_used cohort_flood_cplepa, msize(.5) lwidth(.1) lcolor(black) || ///
       rcap us_diff_af_used us_diff_af_used cohort_flood_cplepa, msize(.5) lwidth(.1) lcolor(black) ///
       legend(off) ///
	   title("Acre-feet withdrawn", size(medsmall)) ///
       ytitle("Difference over pre-treatment period", size(small)) graphregion(color(white) margin(zero)) ///
	   ylabel(-300 -150 0 150) ///
	   yscale(extend) ///
	   xtitle("Cohort", size(small)) ///
	   xscale(extend) /// 
	   name(flood_to_cplepa_afused, replace)
	   
*Plot for acres_irr
twoway rbar lqt_diff_acres_irr med_diff_acres_irr cohort_flood_cplepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rbar med_diff_acres_irr uqt_diff_acres_irr cohort_flood_cplepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rspike lqt_diff_acres_irr ls_diff_acres_irr cohort_flood_cplepa, lwidth(.1) lcolor(black) || ///
       rspike uqt_diff_acres_irr us_diff_acres_irr cohort_flood_cplepa, lwidth(.1) lcolor(black) || ///
       rcap ls_diff_acres_irr ls_diff_acres_irr cohort_flood_cplepa, msize(.5) lwidth(.1) lcolor(black) || ///
       rcap us_diff_acres_irr us_diff_acres_irr cohort_flood_cplepa, msize(.5) lwidth(.1) lcolor(black) ///
       legend(off) ///
	   title("Acres irrigated", size(medsmall)) ///
       ytitle("Difference over pre-treatment period") graphregion(color(white) margin(zero)) ///
	   xtitle("Cohort") ///
	   xscale(extend) /// 
	   ylabel(-200 -150 -100 -50 0 50) ///
	   yscale(extend) ///
	   name(flood_to_cplepa_acresirr, replace)
	   
*Plot for depth_applied
twoway rbar lqt_diff_depth_applied med_diff_depth_applied cohort_flood_cplepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rbar med_diff_depth_applied uqt_diff_depth_applied cohort_flood_cplepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rspike lqt_diff_depth_applied ls_diff_depth_applied cohort_flood_cplepa, lwidth(.1) lcolor(black) || ///
       rspike uqt_diff_depth_applied us_diff_depth_applied cohort_flood_cplepa, lwidth(.1) lcolor(black) || ///
       rcap ls_diff_depth_applied ls_diff_depth_applied cohort_flood_cplepa, msize(.5) lwidth(.1) lcolor(black) || ///
       rcap us_diff_depth_applied us_diff_depth_applied cohort_flood_cplepa, msize(.5) lwidth(.1) lcolor(black) ///
       legend(off) ///
	   title("Depth-applied", size(medsmall)) ///
       ytitle("Difference over pre-treatment period") graphregion(color(white) margin(zero)) ///
	   xtitle("Cohort") ///
	   xscale(extend) /// 
	   name(flood_to_cplepa_depthapplied, replace)

/* Combine the graphs */
grc1leg2 flood_to_cplepa_afused flood_to_cplepa_acresirr flood_to_cplepa_depthapplied, ///
	rows(3) cols(1) loff ///
	xtob1title xtitlefrom(flood_to_cplepa_afused) ///
	ytol1title ytitlefrom(flood_to_cplepa_afused) ///
	labsize(vsmall) graphregion(color(white)) xsize(6.5) ysize(7.5) iscale(.8) ring(1) noautosize
graph export "${dr_output_app}\figureA4.tif", width(6500) height(7700) replace

*********************** Figure A.5 *********************************************
**Center pivot to LEPA
import delimited using "$dr_temp\cp_lepa_prepped.csv", clear
*Drop not-yet-treated during the times when they are treated
drop if wua_year >= cohort_cp_lepa & cohort_cp_lepa!=0
tab treat_status
tab cohort_cp_lepa

*Make mean of each variables
foreach var of varlist af_used acres_irr depth_applied {
egen mean_`var' = mean(`var')
}

*Create long difference variables
xtset wr_group wua_year

*Create max_year variable for the never-treated
bysort wr_group: egen max_year = max(wua_year)
gen max_nt_year = . 
replace max_nt_year = max_year if cohort_cp_lepa==0

*Keep minimum year in dataset and last year before treatment for adopters
	*Keep minimum year and max year for never-treated
keep if wua_year==min_year | wua_year==cohort_cp_lepa - 1 | wua_year == max_nt_year
drop if min_year==cohort_cp_lepa - 1
gen ind_min_year = 0
replace ind_min_year = 1 if min_year==wua_year

foreach var of varlist af_used acres_irr depth_applied {
	gen min_year_`var' = `var'*ind_min_year
	bysort wr_group: egen min_`var' = max(min_year_`var')
	drop min_year_`var'
}

*Keeper years
gen keeper_year = 0
replace keeper_year = 1 if wua_year==cohort_cp_lepa - 1 & cohort_cp_lepa!=0
replace keeper_year = 1 if wua_year == max_nt_year
keep if keeper_year==1

foreach var of varlist af_used acres_irr depth_applied {
	gen diff_`var' = `var' - min_`var'
	replace diff_`var' = 100*diff_`var'/mean_`var'
}

// graph box diff_af_used, over(cohort_flood_cplepa)
// graph box diff_acres_irr, over(cohort_flood_cplepa)
// graph box diff_depth_applied, over(cohort_flood_cplepa)

foreach var of varlist diff_af_used diff_acres_irr diff_depth_applied {
* Use egen to generate the median, quartiles, interquartile range (IQR), and mean. 
bysort cohort_cp_lepa: egen med_`var' = median(`var')
by cohort_cp_lepa: egen lqt_`var' = pctile(`var'), p(25)
by cohort_cp_lepa: egen uqt_`var' = pctile(`var'), p(75)
by cohort_cp_lepa: egen iqr_`var' = iqr(`var')
by cohort_cp_lepa: egen mean_`var' = mean(`var')

* Find the lowest value that is more than lqt - 1.5 iqr 
* this is used to form the lower "whisker" of the boxplot.
gen l_`var' = `var' if(`var' >= lqt_`var'-1.5*iqr_`var')
by cohort_cp_lepa: egen ls_`var' = min(l_`var')

* Find the highest value that is less than uqt + 1.5 iqr 
* this is used to form the upper "whisker" of the boxplot.
gen u_`var' = `var' if(`var' <= uqt_`var'+1.5*iqr_`var')
by cohort_cp_lepa: egen us_`var' = max(u_`var')

* Find any outliers (i.e. values more than 1.5 IQRs from the upper and lower quartiles). 
gen outliers_`var' = `var' if(`var'<=lqt_`var'-1.5*iqr_`var' | `var'>=uqt_`var'+1.5*iqr_`var')
}

replace cohort_cp_lepa = 1990 if cohort_cp_lepa==0


*Plot for af_used
twoway rbar lqt_diff_af_used med_diff_af_used cohort_cp_lepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rbar med_diff_af_used uqt_diff_af_used cohort_cp_lepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rspike lqt_diff_af_used ls_diff_af_used cohort_cp_lepa, lwidth(.1) lcolor(black) || ///
       rspike uqt_diff_af_used us_diff_af_used cohort_cp_lepa, lwidth(.1) lcolor(black) || ///
       rcap ls_diff_af_used ls_diff_af_used cohort_cp_lepa, msize(.5) lwidth(.1) lcolor(black) || ///
       rcap us_diff_af_used us_diff_af_used cohort_cp_lepa, msize(.5) lwidth(.1) lcolor(black) ///
       legend(off) ///
	   title("Acre-feet withdrawn", size(medsmall)) ///
       ytitle("Difference over pre-treatment period (% change)", size(small)) graphregion(color(white) margin(zero)) ///
	   xtitle("Cohort", size(small)) ///
	   xlabel(1990 "Never adopters" 1995 "1995" 2000 "2000" 2005 "2005" 2010 "2010" 2015 "2015") ///
	   xscale(extend) ///
	   name(cp_to_lepa_afused, replace)
	   
*Plot for acres_irr
twoway rbar lqt_diff_acres_irr med_diff_acres_irr cohort_cp_lepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rbar med_diff_acres_irr uqt_diff_acres_irr cohort_cp_lepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rspike lqt_diff_acres_irr ls_diff_acres_irr cohort_cp_lepa, lwidth(.1) lcolor(black) || ///
       rspike uqt_diff_acres_irr us_diff_acres_irr cohort_cp_lepa, lwidth(.1) lcolor(black) || ///
       rcap ls_diff_acres_irr ls_diff_acres_irr cohort_cp_lepa, msize(.5) lwidth(.1) lcolor(black) || ///
       rcap us_diff_acres_irr us_diff_acres_irr cohort_cp_lepa, msize(.5) lwidth(.1) lcolor(black) ///
       legend(off) ///
	   title("Acres irrigated", size(medsmall)) ///
       ytitle("Difference over pre-treatment period") graphregion(color(white) margin(zero)) ///
	   xtitle("Cohort") ///
	   xlabel(1990 "Never adopters" 1995 "1995" 2000 "2000" 2005 "2005" 2010 "2010" 2015 "2015") ///
	   xscale(extend) ///
	   name(cp_to_lepa_acresirr, replace)
	   
*Plot for depth_applied
twoway rbar lqt_diff_depth_applied med_diff_depth_applied cohort_cp_lepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rbar med_diff_depth_applied uqt_diff_depth_applied cohort_cp_lepa, fcolor(gs12) lcolor(black) barw(.5) || ///
       rspike lqt_diff_depth_applied ls_diff_depth_applied cohort_cp_lepa, lwidth(.1) lcolor(black) || ///
       rspike uqt_diff_depth_applied us_diff_depth_applied cohort_cp_lepa, lwidth(.1) lcolor(black) || ///
       rcap ls_diff_depth_applied ls_diff_depth_applied cohort_cp_lepa, msize(.5) lwidth(.1) lcolor(black) || ///
       rcap us_diff_depth_applied us_diff_depth_applied cohort_cp_lepa, msize(.5) lwidth(.1) lcolor(black) ///
       legend(off) ///
	   title("Depth-applied", size(medsmall)) ///
       ytitle("Difference over pre-treatment period") graphregion(color(white) margin(zero)) ///
	   xtitle("Cohort") ///
	   xlabel(1990 "Never adopters" 1995 "1995" 2000 "2000" 2005 "2005" 2010 "2010" 2015 "2015") ///
	   xscale(extend) ///
	   name(cp_to_lepa_depthapplied, replace)

/* Combine the graphs */
grc1leg2 cp_to_lepa_afused cp_to_lepa_acresirr cp_to_lepa_depthapplied, ///
	rows(3) cols(1) loff ///
	xtob1title xtitlefrom(cp_to_lepa_afused) ///
	ytol1title ytitlefrom(cp_to_lepa_afused) ///
	labsize(vsmall) graphregion(color(white)) xsize(6.5) ysize(7.5) iscale(.8) ring(1) noautosize
graph export "${dr_output_app}\figureA5.tif", width(6500) height(7700) replace